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ABSTRACT 

In the presence of self-gravity, we investigate the self-similar dynamics of a relativisti- 
cally hot gas with or without shocks in astrophysical processes of stellar core collapse, 
formation of compact objects, and supernova remnants with central voids. The model 
system is taken to be spherically symmetric and the conservation of specific entropy 
along streamlines is adopted for a relativistic hot gas whose energy-momentum relation 
is expressed approximately by £ = cp with e and p being the energy and momentum 
of a particle and c being the speed of light. In terms of equation of state, this leads 
to a polytropic index 7 = 4/3. The conventional polytropic gas oi P = np^ , where 
P is the thermal pressure, p is the mass density, 7 is the polytropic index, and k is 
a global constant, is included in our theoretical model framework. Two qualitatively 
different solution classes arise according to the values of a simple power-law scaling 
index a, each of which is analyzed separately and systematically. With explicit con- 
ditions, all sonic critical lines appear straight. We obtain new asymptotic solutions 
that exist only for 7 — 4/3. Global and asymptotic solutions in various limits as well 
as eigensolutions across sonic critical lines are derived analytically and numerically 
with or without shocks. By specific entropy conservation along streamlines, we extend 
the analysis of Goldreich & Weber for a distribution of variable specific entropy with 
time t and radius r and discuss consequences in the context of a homologous core 
collapse prior to supernovae. As an alternative rebound shock model, we construct an 
Einstein-de Sitter explosion with shock connections with various outer flows includ- 
ing a static outer part of a singular polytropic sphere (SPS). Under the joint action 
of thermal pressure and self-gravity, we can also construct self-similar solutions with 
central spherical voids with sharp density variations along their edges. 

Key words: hydrodynamics — shock waves — stars: formation — stars: interiors 
— stars: winds, outflows — supernovae: general 



1 INTRODUCTION 

Radiation pressure (e.g., Chandrasekhar 1939, 1960; Rybicki 

6 Lightman 1979), trapped neutrino pressure deep in the 
stellar interior of extremely high nuclear density, relativis- 
tically degenerate materials (e.g., Chandrasekhar 1939), ex- 
tremely hot interior materials of stars, and processes likely 
involved in gamma-ray bursts (GRBs) etc. may be approx- 
imated by an equation of state with a polytropic index of 

7 = 4/3. Statistical physics indicates that the state for a sta- 
tionary radiation field with particles of no rest mass such as 
photons is described by a polytropic relation with an index 
7 = 4/3. It is proven that 7 = 4/3 is also a very good ap- 
proximation for relativistically hot particles with negligible 
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rest mass. Moreover for the stellar structure, Chandrasekhar 
(1939) noted that for an infinitesimal uniform expansion or 
contraction of a gas sphere, it involves precisely a polytropic 
process of an index 7 = 4/3. For a static equilibrium config- 
uration and a presumed P = np^ with a globally constant 
K, the virial theorem indicates that 7 < 4/3 situations are 
unstable and 7 = 4/3 corresponds to a transition from un- 
stable to stable configurations as 7 increases. When 7 = 4/3 
for a static equilibrium configuration, the equilibrium con- 
dition is referred to as the Lane-Emden equation with the 
total enclosed mass M being independent of the system ra- 
dius but dependent upon the value of k. 

On the other hand, based on the conventional poly- 
tropic equation of state P — Kp~' , where k is a global con- 
stant and 7 varies from 1 for an isothermal case to 5/3 for an 
adiabatic process of monatomic gas, astrophysicists explored 
properties of hydrodynamic behaviours in diverse contexts. 
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such as star formation, core formation in molecular clouds 
and supernova explosions etc. For catching the basic physics 
and theoretical simplicity, most works on gravitational stel- 
lar core collapses or stellar explosions were usually carried 
out under the spherical symmetry. Hunter (1962) considered 
the stability of an equilibrium system and the collapse pro- 
cess based on a polytropic hydrodynamics. He demonstrated 
how a dynamical instability during the collapse leads to a 
breakup of the spherically symmetric radial inflow of gas. 
Shu (1977) constructed the isothermal expansion- wave col- 
lapse solution (EWCS) with a weak discontinuity; and this 
self-similar hydrodynamic model has been further developed 
in the past three decades, from the isothermal case (e.g., Shu 
1977; Lou & Shen 2004) to the polytropic case (e.g., Suto & 
Silk 1988; Yahil 1983; Lou & Wang 2006), as weU as to the 
logotropic case (e.g., Mclaughlin & Pudritz 1997). Observa- 
tionally, spectral lines of CS, H2CO and other molecules in 
star forming regions show that the single peak of each molec- 
ular line splits into double peaks with the blue peak brighter 
than the red peak, which has been regarded as a support- 
ing evidence to the Shu model (e.g. Zhou 1992; Walker, 
Narayanan & Boss 1994; Myers et al. 1996). It is generally 
expected that shock waves are also involved in self-similar 
collapse or expansion profiles (e.g., Tsai & Hsu 1995; Shu et 
al. 2002; Bian & Lou 2005). 

Note that all above studies were carried out on either 
the isothermal case or the 7 7^ 4/3 polytropic case. In one 
case, the polytropic case of 7 = 4/3 is treated as a lim- 
iting case (Yahil 1983). In contrast, Goldreich & Weber 
(1980) directly considered homologous core collapse for a 
conventional polytropic gas with 7 = 4/3 by making use of 
the time reversal invariance. They concluded that when the 
pressure decreases by a fraction of no more than about 3% 
from a static polytrope in equilibrium, a homologous core 
collapse would occur in the stellar interior. On the other 
hand, numerical simulation of Bethe et al. (1979) indicated 
a fractional reduction of pressure by 26% in order to initi- 
ate a core collapse for a supernova explosion. Goldreich & 
Weber (1980) tried to reduce this difference by introducing 
an inner core of a progenitor star; Yahil (1983) performed 
his polytropic analysis and treated their result as a limit of 
7^4/3+. 

Meanwhile, specific entropy conservation along stream- 
lines does not necessarily mean a constant specific entropy 
everywhere at all times. A more general distribution would 
be a variable specific entropy in time and radius (e.g. Cheng 
1977, 1978). Fatuzzo et al. (2004) introduced a self-similar 
transformation to formulate a more general problem, which 
involves a scaling index a and another index q. The more 
general equation of state appears to be P oc M''p~'. The 
g = case corresponds to the conventional polytropic gas. 
In fact, according to this more general equation of state, the 
conservation of mass implies the conservation of specific en- 
tropy along streamlines. Nevertheless, Fatuzzo et al. (2004) 
mainly focused on the isothermal cases with nonzero inward 
flow speeds far away in molecular clouds (Shen & Lou 2004). 

Our consideration is on a more general polytropic gas 
with 7 — 4/3 with the speciflc entropy conserved along 
streamlines. By a self-similar transformation, we can ap- 
proach the resulting nonlinear ordinary differential equa- 
tions (ODEs) systematically. Solution properties depend on 
the scaling index a. Given a distribution of variable speciflc 



entropy with time and radius, the result of Goldreich & We- 
ber (1980) can be substantially extended. Meanwhile, many 
counterparts of previously known solutions in the isother- 
mal and conventional polytropic cases can also be derived. 
In particular, several new asymptotic solutions unique to 
7 = 4/3 are also obtained. An important and interesting re- 
sult of our analysis is that under the joint action of thermal 
pressure force and self-gravity, a central spherical void can 
form and evolve in a self-similar manner; this is to be com- 
pared with the central spherical void solution of Fillmore 
& Goldreich (1984b) which considered a collection of colli- 
sionless particles under self-gravity in the Einstein-de Sitter 
expanding universe. 

This paper is structured as follows. Nonlinear adia- 
batic hydrodynamic equations in spherical symmetry and 
self-similar transformation are described in Section 2 for a 
polytropic gas with a polytropic index 7 = 4/3. The exten- 
sions of the classical analysis of Goldreich & Weber (1980) 
are presented in Section 3 and further discussed for a ho- 
mologous stellar core collapse in Section 6.1. We mainly fo- 
cus on cases of g = 2/3 for various solution properties such 
as the requirement on the scaling index a, the property of 
scaling invariance, global analytic solutions, the sonic sin- 
gular surface, the straight sonic critical lines, eigensolutions 
across the sonic critical line, various asymptotic solutions, 
and shock jump conditions in Sections 4. We analyze various 
semi-complete numerical solutions with or without shocks 
and corresponding results in Section 5. Finally, we conclude 
and discuss our main results in Section 6. Three Appendices 
A, B and G are included at the end for technical details of 
derivations and analyses. 



2 BASIC NONLINEAR EQUATIONS AND 
SELF-SIMILAR TRANSFORMATION 

As a theoretical model formulation, dynamical processes 
outlined in introduction are governed by the basic nonlinear 
hydrodynamic equations under the assumption of spherical 
symmetry. We naturally adopt the spherical polar coordi- 
nates (r, 0, (j)) in the analysis. The mass conservation is 



— u— - 

dt dr 



and 



dM 
dr 



(1) 



where M(r, t) is the enclosed mass within radius r at time t, 
the mass density p(r, t) is a function of r and t and u{r, t) is 
the radial flow velocity. The above two relations in equation 
III are equivalent to the mass continuity equation 



dp 1 d . 2 N r, 



dt 



dr 



(2) 



The gas motion is governed by the radial momentum equa- 
tion 



du du 
dt dr 



1 dP 
p dr 



GM 



(3) 



where P{r,t) is the thermal gas pressure and G = 6.67 x 
10"** dyne -cm^ g ^ is the gravitational constant. The Pois- 
son equation relating the mass density p and the grav- 
itational potential ^(r,t) is automatically satisfled with 
d^/dr = GM/r"^. Finally, the conservation of speciflc en- 
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tropy s(r, t) along streamlines is simply 



— +u — 

dt dr 







(4) 



where 7 is the polytropic index. We note that P = fcp^ with 
a constant k is just a particular case satisfying equation 
(|4|. Combining the conservation laws of mass and specific 
entropy, the specific entropy s{r,t) can be an arbitrary func- 
tion s = s{M) of the enclosed mass M{r, t). The entropy is a 
statistical quantity associated with a large number of parti- 
cles, it appears that in this situation, the entropy is frozen in 
particles along streamlines. By this consideration, we have 

da ds dM 

dt " dM~dr ^ ^ ' 

with 

d d d , , 

di^dl^^'d-r 

being the total time derivative along a streamline. 
2.1 Self-Similar Transformation 

In order to solve for self-similar solutions from these nonlin- 
ear partial differential equations, we introduce the following 
self-similar transformation to reduce equations ([l]) — ((4]) to 
nonlinear ODEs, namely 

aix) ^ _ m{x) 



X = At°'r , 
v{x) 



P = 



P = 



p{x) 



where a is an important scaling index parameter and A is 
a dimensional constant coefficient to make the independent 
variable x dimensionless. Here, a{x), m{x), v{x), and p{x) 
are functions of x only and are referred to as the reduced 
density, enclosed mass, velocity and pressure, respectively. 
Now with self-similar transformation (0 , equations (O — (U) 
take the form of 

dm 

(ax + V)—— = (3a 
dx 

dm 
dx 

. .dv 1 dp m 

dx a dx x'^ 

+ log 

, ^ . 1 da ^ dv 
a dx dx 

(Fatuzzo et al. 2004; Wang & Lou 2007). Before proceed- 
ing, we note that these equations are invariant under the 
following time reversal transformation, namely 



)m , 


(8) 


2 

• a , 


(9) 


i)v , 


(10) 


7) , 


(11) 


-) 

x/ 


(12) 



r 



t - 
M 



-t , 
M 



P 



(13) 



Therefore any solution can also depict its inverse process as 
long as this process is reversible (e.g., not involving shocks). 
For example, one solution describing a collapse can be also 
utilized to describe an expansion process. More importantly, 
equation ((S)) implies a division of all cases into three classes 
by whether or not scaling parameter a is greater than, equal 
to or less than —2/3; in general, a is required to be negative. 



This requirement of a negative a is not obvious by equa- 
tions (ffl) — (fT2)) . By the asymptotic solutions (|34]) and (|46t 
at large x derived later, it is necessary to require a < for 
convergent solutions at large x. 



3 HOMOLOGOUS CORE COLLAPSES 

We first analyze the case of a — —2/3 precisely which in- 
cludes the classical analysis of Goldreich & Weber (1980). 
Their model was applied to a stellar core collapse under self- 
gravity prior to the core bouncing in the context of super- 
nova explosions. By equations ((8} and (|9}, we simply have 



[ax + v)x a = 



(14) 



The case of a = everywhere at all time would be a trivial 
solution; for nontrivial solution, the radial fiow velocity is 
thus given by 



V — —ax — 2x/3 



(15) 



and then equation (|11|) requires j — a — 2 leading to 
7 = 4/3 precisely. Here v{x) represents an expansion so- 
lution, or a core collapse solution with the time reversal in- 
variance transformation. Meanwhile, equation (|12|l becomes 
automatically satisfied under this transformation, giving no 
further information or constraint. Taking the derivative of 
equation (|10|) with respect to x, we derive 



_l_d_ 

x^ dx 



a dx 



2 

""+3 



(16) 



Now these equations are not complete yet and a more gen- 
eral description of specific entropy distribution as a function 
of X is allowed. In other words, we already get P oc p*^^ but 
do not know the proportional coefficient as a function of 
(r, t) which is associated with the enclosed mass M{r, i). 
In fact, this point can also been seen by directly compar- 
ing P and p in self-similar transformation (0. Physically, 
log(P/p'*'''') is proportional to the specific entropy s{r,t) in 
a polytropic gas. Once we know the distribution of specific 
entropy as a function of x, the self-similar polytropic flow is 
then determined. 

The flrst cut is to take a constant specific entropy ev- 
erywhere at all times, i.e., P = Kp^^^ with k being a global 
constant. In fact, this is exactly what Goldreich & Weber 
(1980) did. For A = (47rG/K^)^/'' in self-similar transforma- 
tion (O, we immediately obtain p = q*^^ and a second-order 
ODE for a from equation (|16p . We may write a — for the 
convenience of comparison and the second-order nonlinear 
ODE for /(x) with a central condition is then 



dx^ X dx 
/'(0)=0, 



(17) 



where /(O) > is related to the central mass density and 
is an adjustable parameter up to fc. In numerical integra- 
tions, we may encounter f{x) = at a finite a; > under 
certain conditions. If this is the case, an outer travelling 
boundary of the flow system exists. It is fairly straightfor- 
ward to solve equation (|17[l numerically by the standard 
Runge-Kutta scheme (e.g.. Press et al. 1986). 

We have just summarized essential results of Goldreich 
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Figure 1. Profiles of fix), closely related to the reduced mass 
density by a{x) = f^{x), with different values of /(O) are shown 
for a globally constant specific entropy distribution. Here, f{x) is 
normalized by /(O) which is related to the central mass density. 
There exists a limiting value /(O) = fc such that if and only if 
/(O) > fc, the solution curve f{x) vanishes at a finite x > 0. 
The value fc = 4.67047 is numerically determined. To show this 
transition, we take /(O) = 3.0, fc, 5.0, 10.0 in turn as examples 
of illustration. For /(O) < fc, the solution would have a density 
profile of an infinite extent and the radial flow velocity diverges 
for large x. The curves f{x) give self-similar profiles of density 
distribution for a spherical expansion. The process of collapse 
can be also described by a time reversal operation. 



& Weber (1980) in their analysis. Through numerical explo- 
ration, we also find that there exists a limiting value for /(O) 
denoted by fc = 4.67047 (see Figure [ij. With /(O) greater 
than this critical value fc at x = 0, the solution of f{x) 
is confined by a finite x and has plausible physical proper- 
ties. For /(O) < fc, the mass density does not vanish at a 
finite a; > and does not approach zero for large x either. 
By comparing our adjustable parameter /(O) with parame- 
ter A of Goldreich & Weber (1980), we readily establish the 
following simple conversion relation 



A,r(0) = 2/3 



(18) 



Parameter A has a maximum value as noted by Goldreich 
& Weber (1980); their maximum value = 0.00654376 
corresponds to our fc well. Our fc = 4.67047 gives a Am = 
0.00654375. When /(O) is lower than this minimum value 
fc, there is no self-similar solution with vanishing density at 
a finite x 0. Physically, fc corresponds to the minimum 
central density for a homologous or self-similar core collapse 
to be possible. 

In addition to the preceding analysis, our polytropic 
model analysis here does not necessarily require a constant 
specific entropy everywhere (in time and space) and there- 
fore substantially generalizes the work of Goldreich & Weber 
(1980). In fact, we can allow for a fairly arbitrary distribu- 
tion of specific entropy and therefore accommodate a broad 
class of solutions for the density profile. A proper distribu- 
tion of specific entropy can be described by 



where g{x) is a sensible but otherwise arbitrary function. In 
fact, the case studied by Goldreich & Weber (1980) simply 
corresponds to g(x) = 1. For a more general g{x), we readily 
derive a second-order nonlinear ODE for f{x) and a central 
condition, namely 

igf" + (59' + ^) f + -fg' + fa" + = | , 

\ X J X 3 

p'(0) = ^ 3'(0)/(0) + 4g(0)/'(0) = , (20) 

where prime "/" indicates a derivative with respect to x and 
a{x) = f^{x) is the reduced mass density. Now given a value 
of /(O), related to the central mass density, we can solve 
f{x) numerically to determine the self-similar mass density 
profile. The intersection of f{x) with the x axis is the moving 
'boundary' of the fiow system, denoted by xt- 

As we know the density profile, we can calculate the 
enclosed mass m{x) and the ratio between the central and 
mean densities. As shown by equation (|9}, the enclosed mass 
is 

(21) 



(22) 



m{xi,) = / X adx . 
Jo 

Using equation (I16|l . one can readily get 

m{xt) = 2x1/9 ~ 4:xtg{xt)f'{xi) 

and the dimensional enclosed mass is expressed as 



M = 



m(xb) _ 1 
~AHT ^ A^G 



-x'i - 4:xlg{xt)f'{xt) 



The ratio between the mean and central densities is 



P_ 

Pc 



fHO) 



'ig{xt)f'{xi,) 



Xb 



(23) 



(24) 



We are now in a position to make a comparison. With 
g{x) = 1, Goldreich & Weber (1980; GW) computed the 
above quantities; within numerical errors, the results of 
theirs and ours are mutually consistent. Our result of p/ pc 
varies between 0.0066 and 0.0185, while theirs varies be- 
tween 0.0065 and 0.0185. The value of m{xb), similar to 
rlp/pc in GW, increases by a factor of 1.045 when /(O) in- 
creases from fc to a sufficiently large value, which is also 
equal to that of GW. 

As examples of illustration, we shall prescribe specific 
functional forms of g{x) and analyze corresponding solutions 
of f{x) presently in Section 6.3. 



4 VARIOUS SOLUTION PROPERTIES 

In this section, we mainly focus on cases with a 7^ —2/3. 
In these cases, it is still possible for 7 — 4/3 which was not 
considered by Goldreich & Weber (1980) and Yahil (1983). 



4.1 A Preliminary Consideration 

We turn to reduced nonlinear ODEs (|8)) — (|12|) to start our 
discussion. First, a combination of equations ([S} and (O 
immediately gives the reduced mass as 



m{x) 



{ax + v) 2 . s 
-X a(x) 



(25) 



g(x)a'^^^ 



(19) 



(3a + 2) ' 

By equation (|25|l . no confined solution for a{x) by a finite 
value of a; > exists because q = at a finite x directly 
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leads to m = 0, i.e., no enclosed mass at all within this x 
where the mass density vanishes. As a result, no solutions 
can be confined by a finite x > 0. In these cases, the x 
range of both analytical and numerical solutions is infinite 
and some sensible cutoffs need to be introduced for astro- 
physical applications. Moreover, the enclosed mass should 
be always positive such that (ax + v)/{3a + 2) > 0. For 
a < —2/3, we must require v < —ax, while for a > —2/3, 
we should have v > —ax. This is a strict constraint of self- 
similar transformation such that no decreasing solutions of 
v{x) exist for a > —2/3. In general, a should be always less 
than because of the requirement of a real physical system. 
Dividing equation (lllf) by equation ((HJ, we obtain 



p — Com'^a"' 



(26) 



where another index parameter q is defined by 

gEE 2(2 + a-7)/(3a + 2) (27) 

and Co is a constant of integration. This carries an apparent 
physical meaning. It is mentioned earlier that if the specific 
entropy is a function of the enclosed mass, the conservation 
of specific entropy along streamlines is automatically satis- 
fied. The similarity transformation then gives a more specific 
constraint on the form of this function, which is proportional 
to M'' . Substituting the dimensionless quantities for dimen- 
sional ones, we explicitly obtain 



P = CoA^^"^(47rG)^-'G"A/V^ 



(28) 



Here we have two constant coefficients: A is introduced in 
the transformation and Co is a constant of integration. For 
7 7^ 4/3 and thus q 7^ 2/3, we can adjust the value of Co and 
A such that Co = 1 (see Lou & Wang 2007 for more details) . 
In this paper, however, we focus on the case of 7 = 4/3 and 
thus q = 2/3. The constant A no longer plays an important 
role because the exponent index vanishes in expression (|28p 
and thus disappears. In contrast, Co becomes vital in our 
case under consideration. On one hand, the local specific 
entropy is 

s = log (^-^^ = log Co + i log(4^G) + I log(GA/) . (29) 

The value of Co is related to the specific entropy. On the 
other hand, the local polytropic sound speed is 

3p; "V 3 



dP 



1/2 



(30) 



which is also related to the value of Co. The value of Co 
will affect our equations and thus solutions in a nontrivial 
manner. 



Substituting equations (|26|) and pSI) into equations (flOl 
and p2[) . we readily obtain two coupled nonlinear ODEs 



, ^dv 4Co 4/3 /ax 4- da 



ax V 
3a 2 



a -f 



, , 2Co ( ax + v\ 4/3 , , 

2(1-^). (32) 



dv [ax + v) dee 
dx a dx 



Explicit expressions of these two equations for dv/dx and 
da/dx are contained in Appendix |X] Our subsequent anal- 
ysis is based on these two coupled nonlinear ODEs (|31|l and 



Before a further discussion, one notes that besides the 
time reversal invariance, ODEs (|3H) and (|32} are also invari- 
ant under the following scaling transformation, namely 



■ -qx 



3 

V m , 
2 

V P 



(33) 



where i] is an arbitrary positive constant. This scale invari- 
ance only exists when 7 = 4/3 or g = 2/3 and brings us 
considerable convenience in theoretical analysis. 



4.2 Global Analytic Solutions 

Previously, two kinds of analytic solutions were found, 
namely, the static singular polytropic sphere (SPS) solution 
and the Einstein-de Sitter expansion solution in the New- 
tonian regime (e.g., Wang & Lou 2007). We confirm that 
for the current special case of 7 = 4/3, these two solutions 
still exist with certain modifications and constraints. For 
the former, we note that no static SPS solution exists for 
> a > —2/3 because of inequality v > —ax > 0. For 
a < —2/3, we can set n = in the two ODEs and obtain 



Co 



D 2/a 

a — Bx ' 

4/3 



(3a + 2) 
2(a-|-2) V3a + 2 



(34) 
(35) 



where B > is an arbitrary positive coefficient. Unlike pre- 
vious polytropic models with 7 7^ 4/3, parameter Co here is 
specifically determined by a chosen a < —2/3 in the model. 
It implies that the system requires a special relationship 
between the thermal gas pressure and the combination of 
^2/3^4/3 keep the system in a radial force balance. 

The so-called Einstein-de Sitter solution in the Newto- 
nian approximation with a constant mass density also exists 
here for 7 = 4/3. By taking a constant density, we obtain 



2 

v = -x 



-(l + 2v^Co 



^(l + 2^Co)-^ 



(36) 



where Co > is fairly arbitrary. This solution is independent 
of a value as long as a < —2/3 and describes a homogeneous 
expansion in the Newtonian cosmology. In our case, the con- 
stant a is somewhat different from those of the cases with 
g = and 7 7^ 4/3 (i.e., a conventional polytropic gas with 
a constant specific entropy everywhere). We also find that 
this kind of solutions exists only in two situations: one is 
g = and a < —2/3 (see equations (24) and (25) of Fatuzzo 
et al. (2004)), while the other is g = 2/3 also with a < -2/3 
obtained above. 

It is easy to prove that the Einstein-de Sitter solution 
only exists in two cases for g = and g = 2/3. For a being 
constant in ODEs ((8)l — (fT2)) . equation ([9]) gives 



: x^a/3 



(37) 



where a natural boundary condition is simply m(0) = 0. 
Equation ((Sjl then leads to t; = 2x/3 which also satisfies 
equation It follows from ODEs ((STj) and ([lO} that 



Q^Cog(a/3)« 



= (2/9-Q/3)a 



(38) 



with a being a constant. It is clear that for g = 0, we have 
Q — 2/3, while for g — 2/3, we have a different constant 
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a = 2/(3 + 2Co3''''^) as indicated by solution ([Mil at 7 = 4/3 
and a < —2/3. For q not equal to these two special values, no 
Einstein-de Sitter solution is possible. However for a = —2/3 
precisely as in Section 3, equation is of a = form 
and equation (llOp defines a particular form of g{x), that is, 



(39) 



where f{x) = a ' is also a constant. By the property of 
time reversal invariance, this may also describe a particular 
homologous collapse with a{x) being constant and a finite 
reference radius. 



4.3 Singular Surface and Sonic Critical Line 

When the determinant of the coefficient matrix of equations 
l(3T|l and (|32|) vanishes (see Appendix |A]) , i.e.. 



/ , n2 4Co / ax + v\^^^ 4/3 

(ax + v) — X a , 

^ ' 3 V 3a + 2 y 



(40) 



where the right-hand side is the polytropic sound speed 
squared, the relevant ODEs (I3ip and (|32p become singular 
and no finite first derivatives can be obtained (see Appendix 
A). This singularity determines a particular surface, referred 
to as the sonic singular surface in the three variable space 
of X, V and a. Any solution encountering this sonic singu- 
lar surface would diverge except for certain special cases 
which call for additional requirements. One possibility is to 
go across the sonic critical curve with weak discontinuities 
(e.g., Whitworth & Summers 1985 for an isothermal gas) 
or to jump across the sonic singular surface with shocks 
(e.g., Tsai & Hsu 1985; Shu et al. 2002; Bian & Lou 2005; 
Yu, Lou, Bian & Wu 2006; Lou & Gao 2006; Lou & Wang 
2006). Another possibility is to go across the sonic critical 
curve smoothly, for which the values of a and v as well as 
corresponding derivatives satisfies critical conditions at the 
intersection point with the sonic singular surface. 



4.3.1 Determination of the Sonic Critical Line 

A necessary condition for the existence of first derivatives 
a'{x) and v'{x) is to require 



ax + V 
3a + 2' 



(a + l)v- 



2Co 
3 



ax + V 
3a + 2 

2{ax + I 



^1/3 



4/3 

X ' a 



(41) 



This equation defines a unique curve on the sonic singular 
surface, which can be crossed by analytically smooth solu- 
tions; this curve is referred to as the sonic critical curve 
because it is physically related to the local sound speed Cs. 

Reliable numerical experiments can give us valuable 
guidance for conceptual and analytical analysis. Through 
extensive numerical exploration (see Figure it is shown 
that the sonic critical curve defined as the intersection of 
two surfaces given by equations (|40p and (|4ip seems to be 
straight lines starting from the origin in the — i; versus x 
plane and a remains constant along the straight sonic crit- 
ical lines. Using equation (|40p to eliminate a in equation 




Figure 2. Behaviours of the sonic critical lines with different 
values of Co and a. From (a = —0.7, Co = 1.5) above to (a = 
—0.66, Co = 0.5) below in order, the corresponding values of slope 
k are -6.1045, -2.3754, -0.5863, 0.1815, 0.9399, respectively. 
They are all rays starting from the origin in the semi-complete 
solution space of —v versus x. 



(|41l) . a simple derivation gives v 
ing determined by 

7/3 

-f(a + l)A: 



kx with the slope k be- 



3(3a-f 2)^ f a + k 



4Co \3a + 2, 

-(3a + 2)(a + A:)/2 = 2(1 - k){a + k) 



(42) 



Then, the corresponding value of constant a is given by 



3(3a + 2)^ f a + k 



4Co 



3a -I- 2 



4/3 



(43) 



Once k is solved numerically with given values of parameter 
pair a and Co, the relevant sonic critical line is determined 
with a corresponding constant a. As a consistent confirma- 
tion, we can also use equation (|40p to eliminate v and then 
obtain an algebraic equation for a independent of x. The 
same conclusion can be reached. Our extensive numerical 
experiments also agree with our analytical analysis as ex- 
pected (see Figure [2}. 

For isothermal cases, the projection of the sonic singu- 
lar surface coincides with that of the sonic critical line. Thus 
in the —v versus x plane, the behaviour of critical line can 
also show that of the singular surface. Nevertheless, in our 
situation the shape of the singular surface is fairly compli- 
cated and the critical curve is just a special curve embedded 
in it. The projection of the curve in the —v versus x plane 
cannot show the exact shape of the entire singular surface. 
This is very important in the discussion of shocks because a 
shock solution needs to jump across the sonic singular sur- 
face rather than the sonic critical curve. 



4.3.2 Eigensolutions across the Sonic Critical Line 

One solution seldom crosses the sonic singular surface 
smoothly even if it meets the sonic critical line. There are 
also some constraints on derivatives of proper solutions. For 
an arbitrary point along the sonic critical line, denoted by Xc 
here, we can expand an analytic solution in terms of Taylor 
series expansion in the vicinity of this sonic critical point. 
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Because the nonlinear ODEs is of second-order, only the first 
two terms of the series expansion need to be considered. We 
write 



X = Xc + S 



Vc + Svi 



(44) 



where Vc and Oc are the values at the sonic critical point 
with vi and a\ being the corresponding first derivatives of 
V and Q at a;c, and 5 is a small displacement away from the 
critical point Xc- Substituting expression (|44p into coupled 
nonlinear ODEs pf p and (|32p and keeping in mind Vc = kxc 
and expression (|43p for Qc, it is fairly straightforward to 
derive a quadratic equation for vi, namely 



7 2 
3 ' 



7 + 5a - -kj vi - —k + afc + 6fc 

3 2^ „ flc [a + 2(l-fc)] 

— -a — 6a — 6 7- — =0 

2 (3a + 2) 



(45) 



Once the two first derivatives are determined, higher-order 
derivatives can be calculated systematically according to 
nonlinear ODEs (|3ip and p2p . The two eigenvalues of the 
velocity first derivative v\ = dv/dx are independent of the 
position of critical point shown in equation (|45p , which 

is fundamentally related to the scaling invariance equation 
(|33p discussed earlier. Generally speaking, quadratic equa- 
tion (|45p has two real roots or a pair of complex conjugate 
roots. Because of the implicit expression of k, we are not 
able to give an analytical analysis to decide the existence of 
real roots. However in our numerical exploration, all roots 
are real; that is, for a given values of a and Co in our exper- 
iments, the eigenvalue problem has two real roots so far. 

Numerical tests also show qualitatively different be- 
haviours corresponding to the two eigensolutions across the 
sonic critical curve. To examine global properties of these 
two, one can integrate from the critical point both outwards 
and inwards, with initial values given by the series expan- 
sion solutions in the vicinity of the critical point. We call 
the eigensolution which diverges as x approaching infinity 
as Type I while the other that converges at large x is re- 
ferred to as Type II. The classification of Types I and II 
solutions is merely for the convenience of discussion. 



4.4 Various Kinds of Asymptotic Solutions 

In addition to the two global analytic solutions presented 
above, we also derive various asymptotic solutions either 
near the centre (i.e., x — > 0"*") or at infinity (i.e., x —> -l-oo). 
Because a > —2/3 strongly limits solution behaviours such 
that V remains always positive and diverges as x approaches 
infinity, we focus on cases of a < —2/3. Previously known 
asymptotic solutions will guide us to search for their counter- 
parts in the case of 7 = 4/3 and other possible new solutions 
should they exist. 

By assuming \v{x)\ and a{x) to be nonincreasing func- 
tions for large x, the first typical asymptotic solutions are 



Hx 



2/a 



Lx 



(a+l)/a 



+ Kx 



(2+a)/a 



(X > 1) (46) 



where H > and L are two constants of integration, 7 = 
4/3, and K is determined by 



K = 



aH 



(3a + 2) 



2CoH 



3a 



-1/3 



(a 4- 2) 
(3a + 2) 



(47) 



The free parameter L was first obtained by Whitworth & 
Summers (1985) for an isothermal gas fiow. Cases of L = 
correspond to asymptotic breeze {K > 0) or contraction 
{K < 0) solutions, depending on whether 11(3:) is positive 
or negative at large x. For the first leading term of v{x) in 
dimensional flow velocity 



(48) 



it gives a background fiow at infinity and a convergent fiow 
speed should be required such that — 1 < a < 0. The sign 
of L decides the asymptotic flow direction. Cases of L > 
correspond to outflow or wind solutions while cases of L < 
correspond to contraction or inflow solutions (Lou & Shen 
2004; Lou & Wang 2006, 2007). 

In addition, another asymptotic solution at large x is de- 
scribed below; this asymptotic solution may be regarded as 
a perturbation to the exact Einstein-de Sitter solution (|36p . 
Assuming a series expansion solution approaching Einstein- 
de Sitter solution (|36p as a; +00, we write 



V = -X + Ex 



0+1 



-2v^Co 



+ 0{x''+') , 

+ Fxi^ + o(^x'^^ , 



(49) 
(50) 



where E and F are two constants to be determined and the (3 
parameter is required to be negative, i.e., /3 < 0. With x 
+00, nonlinear ODEs pip and (|32} lead to the following 
linear homogeneous algebraic equations for E and F, namely 



{P + 3)aoE+ (a+ - ](5F = , (51) 



"+0^+(3^r2) 



^CoP+l + ^-^Co]F, (52) 



where ao = (2/3)/(l -I- 2y^Co) is the Einstein-de Sitter 
constant reduced density. For nontrivial solutions of E and 
F, the determinant of the coefflcients in linear equations 
(|5ip and (|52p for E and F should vanish. By this condition, 
a quadratic equation of /3 appears, namely 

[(3a + 2)2 -4^Coao]/?^ 

20v^Coao)/3- 6 = . 



+ 3a + 2 



(53) 



After solving quadratic equation (|53p . we retain the negative 
roots of P for the consistency of our approximation. We shall 
show numerical solution examples that approach asymptotic 
solutions (|49p and (|50p at large x presently. 

We have just examined asymptotic behaviours of solu- 
tions for large x. We now turn to asymptotic solution be- 
haviours in the regime of small x. First, we flnd a core col- 
lapse solution as the counterpart of the isothermal free-fall 
solution of Shu (1977), namely 



m(0) 



-,1/2 



j(3a + 2)| 



-[2m(0)]i/2 
[2m(0)]i 



-3/2 



for a < -2/3 , 
for > a > -2/3 



(54) 



where m(0) is the limit of 771(2;) as x goes to 0^. Because 
both mass density and fiow velocity are divergent as x O""" , 
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a physical cutoff needs to be set somewhere at a small x 
as the inner reference 'boundary' of the flow system under 
consideration, or regarded as a reference surface surrounding 
a central compact object. Comparing gravity and thermal 
pressure force, we immediately have 

1 dP m 1 dp 
p dr a dx 



GM 



-l-3{7-l)/2 



(55) 



for a; ^ 0^. Near the centre, the gravity becomes very much 
stronger than the thermal pressure force and therefore dom- 
inates in the process of gravitational core collapse. Accret- 
ing materials then fall towards the centre in the form of 
an almost free fall unimpeded by pressure. This asymptotic 
solution represents such a physical scenario that materials 
accelerate to fall towards the centre under the overwhelm- 
ing self-gravity so that particles gain increasing speed and 
acceleration to impact the central object. For black holes, 
accreting materials are absorbed more effectively. 

As to the so-called Larson-Penston (LP) type solutions 
at small x (Larson 1969a, b; Penston 1969) with no flow and 
flnite density at the centre, we can show that the existence 
of LP type solutions of 7 = 4/3 requires special conditions 
(Lou & Shi 2007 in preparation). With a LP solution in the 
form of a Taylor series expansion near the centre, namely 



v(x) = '^v„x", a(x) = I 



(56) 



where index n runs through non-negative integers, nonlinear 
ODEs (|3ip and (1321) require the constant term vo to be zero. 
After straightforward calculations of ODEs by substituting 
v{x) and a(x), we just obtain Einstein-de Sitter solution 
(|36p . Thus solutions with both flnite density and velocity 
(including the LP type solutions) near the centre may only 
exist under rare situations (see Appendix [B] for details). 

However, if the mass density diverges instead of being 
flnite at small x, a new asymptotic solution can be derived. 
Let us consider the leading order terms of such solutions in 
the form of 



Rx 



Nx^ 



with A < 



(57) 



where R, N and A are three parameters to be determined. 
Then nonlinear ODEs ((311 and give 



A = 



a + R ' 
2Co(2 + a^2i^) 
a + R 



a + R 
3a + 2 



1/3 



= 



(58) 
(59) 



with > being arbitrary. The Type II eigensolution men- 
tioned above just approaches this kind of new asymptotic 
solution at small x. For a < —2/3, we require R < —a in 
order to keep the enclosed mass positive and that A < also 
requires R < 2/3. Simple analysis of equation ((59} shows 
that Co has a critical value 2'^^^ /6 ~ 0.4200 (see Appendix 
C for details), below which R has no real root smaller than 
2/3. For Co > 0.4200, there two real roots of R for such 
asymptotic solutions at small x. For Co < 0.4200, this kind 
of asymptotic solution does not exist. A possible inference is 
that Type II eigensolutions may be truncated before reach- 
ing the origin x = 0. Later numerical solutions conflrm this 
point. 

This solution represents a situation in which the ther- 



2 da 
X a—— ~ 1 



(60) 



mal pressure force is comparable to the self-gravity, that is 
GM 1 dP _ m 1 dp 

p dr x'^ a dx dx 

As the velocity magnitude decreases at small x, the thermal 
pressure force actually becomes somewhat larger than the 
self-gravity. 

Lou & Wang (2006) found a novel "quasi-static" solu- 
tion for a conventional polytropic gas with 7 7^ 4/3 and 
proposed a rebound shock model for supermova explosions 
(see also Lou & Wang 2007). We flnd the counterpart of this 
"quasi-static" asymptotic solution in our case of 7 = 4/3. 
Static SPS solutions are described by equation (|34|l and we 
then introduce next-order perturbations such that 

v~Vx^ , (61) 

a ~ Bx^'" + Wx" , (62) 

where ^ and a are two exponents and V and W are two 
coefflcients to be determined. Note that parameter Co has 
been specifled by equation (|35|l when discussing static SPS 
solutions. It is natural to require ^ > 1 and a > 2/ a in 
reference to SPS solutions. Substituting expressions (|6H) and 
((52} into nonlinear ODEs (|3IJ and (15^ with a sufficiently 
small X and Co expression (|35|) . we obtain two equations for 
coefficients V and W , namely 



e-l=a-- 
a 

2 \ W 
i+-+2\v + {aG-2)— =0 



2\ W 
2+ - ]V +{aa-2) — 
a J Jd 







(63) 
(64) 

(65) 



For nontrivial solutions of V and W, we must require the 
determinant of equations (|64|l and (|65|) to vanish. The result- 
ing equation together with condition (|63|) lead to a quadratic 
equation of ^ (Lou & Wang 2006). The relevant root of ^ is 



C = -4(l + l/a) 



(66) 



while the other ^ — 2/a root of the quadratic equation is 
unacceptable. As ^ > 1 is required, we then have inequal- 
ity —4/5 < a < —2/3. This appears somewhat different in 
certain aspects of Lou & Wang (2006): (i) index ^ is always 
real (no possibility for a pair of complex conjugate roots) 
and only one root is valid for 7 = 4/3; (ii) occasionally, 
both real roots of this index in Lou & Wang (2006) may 
be valid for 7 7^ 4/3; (iii) this index may become a pair of 
complex conjugate roots for 7 7^ 4/3, leading to asymptotic 
oscillations; and (iv) the allowed range of a = —n is larger 
here for 7 — 4/3. 

4.5 Shock Jump Conditions 

When a faster flow catches up to a slower one, a shock 
wave can form and propagate in stellar winds, molecular 
clouds and stellar interiors. A shock occupies a narrow re- 
gion with discontinuities in density, pressure, temperature, 
entropy and flow velocity (e.g., Landau & Lifshitz 1960). 
In our model framework, we are interested in self-similar 
shocks which are "flxed" in a self-similar proflle (e.g., Sedov 
1959). Besides crossing the sonic singular surface smoothly, 
a flow solution can also jump across it by shocks which ex- 
tends physical solutions with various possibilities. In fact. 
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shock phenomena are ubiquitous in astrophysical systems. 
Across a shock front and in the shock framework of refer- 
ence, conservations of mass, momentum and energy hold, 
and so does the second law of thermodynamics, the increase 
of entropy from the upstream side to the downstream side. 
In this subsection, subscripts 1 and 2 always represent phys- 
ical quantities of upstream and downstream sides of a shock, 
respectively. 

In the shock framework of reference, the three conser- 
vation laws of mass, momentum and energy correspond to 
the following three equations, namely 



Pi(mi 



Pl+pl{ui—Us) = P2 + P2{U2 ~ Us) , 
(Ul — Ms)^/2 + TOl = (lt2 — Ms)^/2 + W2 



(67) 
(68) 
(69) 



where Us is the shock speed in the laboratory framework 
and w denotes the heat function defined by 

7P _ 4P 

~ (7 - 1)P P 



(70) 



for a polytropic gas with 7 = 4/3. We usually introduce 
self-similar transformations separately for the upstream and 
downstream sides of a shock. However unlike previous work 
of Lou & Wang (2006), the parameter A can be arbitrary for 
the case of 7 = 4/3 and g = 2/3 so that the same self-similar 
transformation is valid on both sides of a shock0 Because 
shocks are "fixed" in a self-similar profile, the scaling pa- 
rameter a is unchanged across a shock. After this self-similar 
transformation, the three conservation equations (|67|) 
and (|69[) can be reduced to 



02172 



— + ai6'i = — + a2i 



4pi 



4^ 



X g 
4p2 

a2x'i 



el 



(71) 
(72) 

(73) 



where Xs is the shock location (also related to the shock 
speed), and 6i is defined by Vi/xs -\- a. We now manage to 
solve the downstream parameters from the upstream param- 
eters. Using equations (|71|l and (|72|l to eliminate P2 and 02, 
we obtain a quadratic equation for 62- Because of the in- 
variance for exchanging subscripts 1 and 2, we neglect the 
trivial solution 61 — 62 and obtain 



8pi 



4- ■ 



7ai6iXs 7 



(74) 



where Xs is a chosen value for 7 = 4/3 with no difference 
between the upstream and downstream sides of a shockQ It 
is then straightforward to solve for other downstream quan- 
tities 

a2 = aiei/e2 , (75) 
P2 = Pi + {ai9i - a29l)xl . (76) 

We now know p2 and Q2 and m,{x) is continuous across a 



shock, the downstream coefficient C02 = P2 1 (jn^ ^ ^ p^ '^) is 
then determined. 

Besides the three conservation laws, the second law of 
thermodynamics will check whether this solution is physi- 
cally appropriate. Here the second law is satisfied as long as 
the downstream coefficient C02 is greater than the upstream 
coefficient Coi. It is also convenient to introduce the Mach 
number here 



Mi 



(77) 



where Ci (i — 1, 2) is the sound speed (upstream, down- 
stream). Hence equation (|74p can be rewritten as 



U2 



02 

01 



7Ml ~^ 7 



which is equivalent to 
2 + (7-l)A1? 



Ml 



Ml 



2'yMl - (7 - 1) 8Ml 



(78) 



(79) 



in terms of Mach numbers. The increase of specific entropy 
requires that the pressure of the upstream side is lower than 
that of the downstream side; this leads to several inequalities 
below 



Ml 



> 1, 



Mi < 1, 



Cl < C2 



(80) 



Qualitatively speaking, the upstream flow is supersonic 
while the downstream flow is subsonic. 

Reciprocally, we can also calculate quantities of the up- 
stream side from those of the downstream side following the 
same derivation procedure for shock conditions (|7H) — (|73p . 
One should note that from equation (|79|l . the physical con- 
straint on A^i is 1/8 < A^2 < 1 for 7 = 4/3. Therefore when 
we calculate upstream variables from downstream variables, 
the self-similar shock position should be chosen within a 
certain sensible range such that Ml falls within this speci- 
fied range. Otherwise, no physical solutions for a self-similar 
shock can be constructed, i.e.. Ml becomes less than unity. 

During the construction of numerical shock solutions, 
we sometimes specify asymptotic solutions at large x and 
integrate inward. In this case, we specify physical variables 
on the upstream side of a shock first and then derive physi- 
cal variables on the downstream side of a shock. We do not 
encounter troubles in choosing shock positions. However, in 
various occasions, we may need to specify asymptotic solu- 
tions at small x and integrate outward. In such situations, 
we specify physical variables on the downstream side of a 
shock first and then derive physical variables on the up- 
stream side of a shock. For these cases, we need to make sure 
that the downstream Mach number satisfies the inequality 
1/8 < Ml < 1 for our chosen shock positions. Otherwise, 
upstream variables may become unphysical as this happens 
in our numerical exploration. 



5 NUMERICAL SOLUTIONS AND RESULTS 



^ In the analysis of Lou & Wang (2006), the coefficient A is re- 
lated to the local sound speeds, which are different for the up- 
stream and downstream sides of a shock. 

^ For 7 ^ 4/3, Xs is generally different on the two sides of a shock 
(e.g., Wang & Lou 2007). 



We have derived global analytical solutions, various asymp- 
totic solutions at both large and small x, two eigensolutions 
across straight sonic critical lines for 7 = 4/3, and self- 
similar shock conditions. We are now in a position to con- 
struct various global semi-complete solutions numerically by 
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Figure 3. Several numerical breeze solutions with L = are 
shown as the solid curves with parameters H = I, 0.5, 0.1, 0.01 
and corresponding m{0) = 17.75, 8.75, 1.70, 0.16 at a = -0.68 
and Co = 0.4. The straight dashed line passing through the origin 
is the sonic critical line. Near the core, all solutions approach the 
asymptotic free-fall behaviour (see equation 154 l l at small x. 

utilizing and matching these solutions. In reference to the 
sonic singular surface, we divide all numerical solutions into 
three classes: those avoiding the sonic singular surface, those 
crossing the straight sonic critical line smoothly and those 
with shocks. We construct and discuss these solutions in or- 
der. 

5.1 Solutions not Crossing the Singular Surface 

It is straightforward to construct global numerical solutions 
without encountering the sonic singular surface. Starting 
with the convergent asymptotic solution (|46p at large x (e.g., 
X = 100 in our numerical experiments), we integrate back 
towards the centre. This procedure works fine unless the 
solution runs into the sonic singular surface. The solutions 
diverge near the centre, approaching the free-fall asymptotic 
solution (|54p as a; — > 0"*". Numerical integrations outwards 
from the centre tend to be unstable in the sense that the 
determination of the two parameters H and L is fairly sen- 
sitive to the value of m(0). The reason is that there is only 
one parameter m(0) to be decided in the inner part while 
the outer part involves two parameters H and L. In numer- 
ical procedures, using two parameters (e.g., H and L in this 
case) to decide one parameter (e.g., m(0) in this case) is 
stable, while an outward integration from small x to large 
X, using one parameter to decide two parameters, tends to 
be sensitive to the numerical accuracy. 

We adjust the solution parameters H and L and the 
relevant parameters a and Co to explore various solutions 
(see Fig. [3]). When a solution goes back towards the origin, 
it matches with asymptotic free-fall solution (|54|) and gives 
the corresponding m(0) value. For a = —1 in expression 
H46p , solutions with L = have vanishing velocities at infin- 
ity; solutions with L > and L < correspond to constant 
outflows and inflows at infinity, respectively. Such solutions 
offer the following scenario: at the beginning time {t = 0), 
the gas system is stationary or has a velocity outwards or 
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Figure 4. A comparison of —v{x) solution with the sound speed 
c{x) in the flow. Relevant parameters are: a = —0.68, Co = 0.478, 
H = 0.5, m(0) = 8.72. The dashed line is the sonic critical line. 
The solid line is the numerical solution of —v{x) while the dash- 
dotted line represents the sound speed in the laboratory frame- 
work. The solution —v(x) goes from subsonic at large x to super- 
sonic at small x. 



Table 1. Parameters adopted for examples of numerical solutions 
are summarized in this Table 1. All these solutions do not en- 
counter the sonic critical line. Using the standard Runge-Kutta 
scheme of fourth order with initial values calculated from the 
asymptotic solutions 1 146 l l at a sufficiently large x (e.g., at x = 100 
in our numerical integrations). Integrating back towards the ori- 
gin, we match the free-fall asymptotic solution l)54|l as a; — » 0+ 
and estimate the values of m(0) relevant to central mass and mass 
accretion rate. 



a 


Co 


H 


K 


L 


m(0) 


-0.68 


0.4 


2.0 


-13.46 





36.0 


-0.68 


0.5 


2.0 


-8.33 





35.8 


-0.68 


0.6 


2.0 


-3.20 





35.6 


-0.68 


0.4 


1.5 


-10.99 





26.8 


-0.68 


0.4 


1.0 


-6.73 





17.7 


-0.68 


0.4 


0.6 


-4.04 





10.5 


-0.68 


0.4 


0.2 


-1.35 





3.44 


-0.68 


0.4 


0.01 


-0.067 





0.162 


-0.67 


0.4 


0.001 


-0.041 





0.0661 


-0.70 


0.4 


1.0 


-1.56 





7.32 


-0.72 


0.4 


1.0 


-0.62 





4.59 


-0.74 


0.4 


1.0 


-0.30 





3.29 


-0.67 


0.4 


0.001 


-0.041 


0.2 


0.0651 


-0.68 


0.5 


1.0 


-4.16 


0.2 


17.6 



inwards with its mass density profile proportional to r^''". 
Under the joint action of self-gravity and thermal pressure 
force, the entire system evolves into a central collapse even- 
tually. Around the central region, the inward self-gravity is 
always larger than the thermal pressure force so that ma- 
terials are accelerated towards the centre. Nothing singular 
happens as the local flow speed reaches the local sound speed 
(see Fig|4|. When approaching the centre, the fluid is almost 
in a free-fall state. Because of our presumed spherical sym- 
metry, something must happen around the centre to destroy 
the similarity flow or spherical symmetry. For example, a 
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Figure 5. The polytropic EWCS with 7 = 4/3. The point of 
weak discontinuity is Xe = 1. The outer part of the system is 
static SPS while the inner part approaches the free-fall asymptotic 
solution. The enclosed mass at the centre m{0) is around 1.30 
and the collapsing mass m(l) is around 1.35, indicating that the 
majority of materials, around 96.3% concentrates in the centre. 
Such polytropic EWCS of 7 = 4/3 exists only for k = with a 
pair of parameters a = —0.68 and Co = 0.6623. 



strong radiation shock may emerge surrounding the centre. 
Or, a black hole may take all accreting materials in. 

The isothermal expansion-wave collapse solution 
(EWCS) of Shu (1977) was regarded as a limit for a fam- 
ily of solutions without encountering the sonic critical line. 
In other words, at a particular critical point along the 
sonic critical line, one of the two eigensolutions leads to the 
static singular isothermal sphere (SIS) while the other leads 
to a solution matching the central free-fall solution. Thus a 
semi-complete global solution with a weak discontinuity is 
constructed, connecting the two eigensolutions at Xe, with a 
static SIS for x > Xe- Similarly, for a particular pair of a and 
Co values (see equation I35|) . a static SPS with 7 = 4/3 ex- 
ists so that we can also construct the counterpart of isother- 
mal EWCS. As every point x is equivalent in the sense of 
the scaling invariance (|33|l . we simply take = 1 with- 
out loss of generality. For this special pair of a = —0.68 
and Co = 0.6623, we have k — for the slope of the sonic 
critical line, and the two corresponding eigensolutions are 
vi ^ and vi = 3(1 + 5a/7) = 1.54. Using wi = with its 
corresponding ai to integrate outwards, we obtain the outer 
part of a SPS as the static outer envelope. Meanwhile, using 
Vi = 3(l-|-5a/7) = 1.54 to integrate back towards the centre, 
we obtain a central free-fall solution. Together, we have con- 
structed a polytropic EWCS with 7 = 4/3 (see Fig. (5)). Let 
us consider the enclosed mass m{x), where m(0) is the point 
mass at the centre and m(l) — m(0) is the mass collapsing 
towards the centre. Our numerical result is m(0) = 1.30 and 
m(l) = 1.35. That is, 96.3% of the total mass concentrates 
in the central object and only 3.7% is collapsing towards the 
centre. 

We can also construct other semi-complete global solu- 
tions without encountering the sonic singular surface. Start- 
ing from quasi-static solution (|61|l and (|62p at small x with 
fc = and V > in equation (|61|l , straightforward numerical 
integrations lead to global semi-complete solutions without 
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Figure 6. Four quasi-static numerical solutions with fc = 
and y > in quasi-static asymptotic solution II6II I for small 
X. The scaling parameter is a = —0.68 and the corresponding 
Co = 0.6623 is computed from equation 13511 . The parameter B 
in the static solution is set to 1. By these numerical solutions, it is 
clear that the larger the V is, the more rapidly the solution is in 
the aberrance of the static state and approaches the Einstein de 
Sitter solution 1491 1 and 1501 1. The straight dashed line is = 2x/3 
for the Einstein-de Sitter expansion. 

encountering the sonic singular surface. Figure [6] shows such 
examples at a = —0.68 with a corresponding Co = 0.6623. 
These solutions never vibrate towards small x according to 
our analysis (i.e., no complex conjugate roots are possible). 
With increasing x, they approach asymptotic expansion so- 
lution (|49p and H50|l rapidly. Numerical results indicate that 
with a V > 0, the gas begins to flow outwards and ap- 
proaches a constant density. The larger the value of V is, 
the stronger the perturbation is, and the more rapidly the 
solution approaches the Einstein-de Sitter expansion phase 
with V — 2x/3. 

Using new asymptotic solution (|57[) at small x, we can 
also construct global semi-complete solutions. As the sonic 
critical line is not enough to describe the relative position 
between numerical solutions and the sonic singular surface, 
we define a velocity Vc such that 

vc = - (7Com^/^Q'''^) - ax , (81) 

to represent a Vc curve on the singular surface, depending 
upon the solution for the reduced enclosed mass 771(2;) and 
mass density a{x) together with adopted parameters Co, 
a and 7 = 4/3. The purpose is to compare solution v{x) 
against Vc for the possibility of encountering the sonic sin- 
gular surface. From this definition, it is easy to see that if a 
solution v{x) meets the sonic singular surface at some point, 
this point must the intersection of the solution curve v{x) 
and the Vc curve thus defined. In equation (f57)) , R, N and A 
are three parameters to be determined. For an appropriate 
combination of these parameter values, the solution may not 
run into the sonic singular surface and eventually converge 
to asymptotic solution (|49|l and (|50p at large x. This solution 
gives the following scenario: the central region is occupied 
by a high-density core with an outside medium of constant 
density moving outward. Once the inner core begins to move 
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Figure 7. An example for the case of a = —0.68 and Co = 0.9 
is illustrated. The corresponding slope k = —0.3745 for the sonic 
critical line. It has two real roots of R: one \s R = —1.1235 and 
the other is R = 0.6535. The inner part can be described by 
asymptotic solution JSTj at small x with other parameters N = 
0.500 and A = —2.9778, while the outer part can be described by 
asymptotic solution II49I I and JSOj at large x with parameter /3 = 
—0.7362. The two dash-dotted curves in the left panel are Vc{x) 
as defined by expression II81I I related to the relevant solutions. 

towards the centre, the outer part decelerates to and also 
begins to move inwards. We also see that the mass density 
decreases rapidly as x becomes large. Figure [7] shows a pair 
of such examples. 

5.2 Solutions crossing the Straight Sonic Critical 
Line Smoothly 

Only eigensolutions along the sonic critical line, derived 
in subsection 4-3-2, can cross the sonic singular surface 
smoothly. To investigate their properties, we start numeri- 
cal integration from the vicinity of a sonic critical point with 
initial values given by one of the relevant eigensolutions. The 
analysis here becomes much simpler in light of the scaling in- 
variance transformation (|33|) and properties for every sonic 
critical point are the same and hence an arbitrary point x 
can be chosen for a certain pair of a and Co parameters. 
By solving for eigensolutions along the sonic critical line 
and then extending the eigensolutions globally by numeri- 
cal integrations, we find that the two types of eigensolutions 
have qualitatively different properties. A type I solution ap- 
proaches the free-fall asymptotic solution (|54|l in the inner 
part for small x and has an outer asymptotic solution de- 
scribed by solution (|49|) and (|50|l at large x. The physical 
scenario is that initially the gas with a constant density has 
a tendency to move outwards, because of the thermal pres- 
sure against gravity. The gravity force competes with the 
thermal pressure and wins eventually as time goes on, and 
then the gas begins to decelerate and accelerate to collapse 
towards the centre. Now the self-gravity dominates the ther- 
mal pressure force completely so that materials approach a 
free fall and finally smash onto the central object. 

In contrast, a type II solution has a quite different be- 
haviour. In the vicinity of the origin, the velocity vanishes 




X 



Figure 8. Two examples of semi-complete solutions crossing the 
straight sonic critical line in two possible eigendirections with pa- 
rameters a = —0.68 and Co = 1.5. The sonic critical line is rep- 
resented by the straight dash-dotted line with a negative slope 
k = —2.3754. The solid curves are the possible eigensolutions 
crossing the sonic critical line, one of which is at i: = 0.5 and the 
other is at X = 1.25. The two eigensolutions approach different in- 
ner and outer asymptotic solutions described previously. Because 
of the scaling invariance of the ODEs, these two eigenproblems 
and the corresponding eigensolutions are actually the same. Take 
the eigenproblem at x = 1.25 for example. Type I solution ap- 
proaches equation 1154 l l with m(0) = 38.4 near the centre and 
has an outer asymptotic solution described by equation (1501 and 
ggjl with parameters /3 = -1.6204, E = 0.086 and F = -0.785. 
Type II solution approaches equation Il57|l at small x with pa- 
rameters A = -3.00, R = -7.90 and N = 0.179. For large x, its 
behaviour can be represented by equation Il46p with H = 2.055 
and L = —11.392. Please note that for a > — 1, —v{x) of a Type 
II solution always vanishes at large x\ in this case of a = —0.68, 
the leading term of —v{x) at large x scales as a:"'^-'*^. 



while the mass density diverges as described by asymptotic 
solution (|57p around small x, while fiow behaviour at large 
X can be described by asymptotic solution (|46|l . Based on 
the value of Co compared with the critical value of 2'^^'^ /Q, 
which determines whether R has real roots, Type II so- 
lutions can be divided into two subtypes. Subtype I: For 
Co > 2*''^/6 corresponding to the existence of real roots 
of R, Type II solutions will approach a; = as described 
by asymptotic solution (|57[) . A special case is the SPS so- 
lution with 7 = 4/3 when k = 0. One can prove that the 
value of Co for A: = is not smaller than 2'^^^ As dis- 
cussed earlier, EWCS with 7 = 4/3 can be constructed here 
by connecting two branches of the two eigensolutions along 
the sonic critical line with = 0. These subtype solutions 
describe the following scenario. The outer part of the fluid 
system has a common flow behaviours which can be an in- 
flow, or an outflow, or even a static envelope, while in the 
inner part, the pressure force and self-gravity compete with 
each other such that the magnitude of the radial speed re- 
mains flnite and eventually vanishes at the centre. Subtype 
II: Co < 2*''^/6 so that asymptotic solution l|57p does not ex- 
ist. A numerical integration backwards would be truncated 
before x becomes sufficiently small. Physically, the enclosed 
mass m(x) is related to the factor aa; -1- w by equation (|25|) . 
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Figure 9. When Co is below 2^/^/6 = 0.42, the asymptotic solu- 
tion l|57|l no longer exists. The type II eigensolution is truncated 
before x = 0. This figure shows an example of such behaviour. 
The relevant parameters are a = —0.68 and Co = 0.3 with the 
corresponding slope k = 0.3571 for the sonic critical line. We 
obtain eigcnsolutions across the sonic critical line at Xc = 1.25. 
Apparently, the enclosed mass becomes before x reaches the 
origin, i.e., in the —v versus x diagram, the velocity curve meets 
the line of ax + v = where the solution is truncated. 



Thus if v{x) curve occasionally approaches the straight line 
ax + V = 0, then there is no material inside this 'radius' 
Xv In other words, a spherical void surrounds the centre 
and expands as time goes on in a self-similar manner. For 
a < —1 cases, the boundary of such a spherical void ex- 
pands with deceleration and the edge radius is proportional 
to t~"-~^. From ODEs the behaviours of density, 

velocity and pressure can be deduced. The enclosed mass 
within that point is zero, the reduced density a is finite, the 
reduced pressure p approaches zero there, and the pressure 
gradient dp/dx = —a{a + l)a remains finite according to 
equation (|10[) . Equation requires the following limit 

{ax + v) da _ 2(2 + a - 7) 



lim 

ax-\-v — *0 



dx 7 
and it follows from equation (I12|l that 

2(2 + a- 7) 



J = 2(l + a 
dx 



7 



(82) 



(83) 



remains finite there. Numerical results also confirm the sit- 
uation that the enclosed mass 771(2;) becomes before x 
reaches the origin. An example of a = —0.68 and Co = 0.3 is 
shown in Figure [9] for the reduced quantities, such as 771(2;), 
a{x) and —v{x) in top, middle, and bottom panels respec- 
tively. This shows the real possibility of a spherical void 
occupying the central region of a certain astrophysical sys- 
tem (e.g., clouds, bubbles, planetary nebulae, stars or super- 
nova remnants etc.) during its evolution under joint action 
of thermal pressure and self-gravity. Previously, Goldreich & 
Fillmore (1984b) discussed coUisionless particles with self- 
gravity in an Einstein-de Sitter expanding universe. Steep 
perturbations can give rise to voids surrounded by overdense 
shells with sharp edges. Our preliminary results here show 
that in addition to the expansion of universe, a spherical 
matter system with thermal pressure against self-gravity can 



Figure 10. Two shock flow solutions are illustrated here. These 
solutions connect free-fall asymptotic solution at the inner part 
(small x) with an inflow or contraction at inflnity (large x). The 
solid curves in the flgure presents solutions and the dashed curves 
are corresponding curves of Vc deflned in equation JSTJ. The scal- 
ing index a = —0.68. The two solutions share the same down- 
stream branch with the free-fall solution parameter m{0) = 0.030 
and C2 = 1.00 and it crosses the sonic critical line smoothly 
at Xc = 0.2. The upstream side shows an inflow described by 
equation II46I I. with a set of parameters {Ci, Xs, H, L}. Two ex- 
ample solutions correspond to {0.5973, 0.40, 0.0019, -0.1322} and 
{0.9580, 0.30, 0.0019, -0.1246}. 



also lead to the formation of a central spherical void with 
an overdense shell along a sharp edge. 



5.3 Self-Similar Flow Solutions with Shocks 

Global behaviours of eigensolutions crossing the sonic crit- 
ical line have been explored numerically. Starting from the 
two eigensolutions on the sonic critical line and integrating 
towards small x, one will approach the free-fall asymptotic 
solution (|54|l and the other will approach the new asymp- 
totic solution (|57|) (see solution examples in Fig. [SJ . Type II 
solutions in Fig. |8] touch the sonic critical line twice. Other 
than this special situation, due to the scaling invariance 
property, we are unable to construct any global solutions 
across the sonic critical line twice smoothly which are pos- 
sible in the isothermal cases of Lou & Shen (2004) and the 
conventional polytropic cases of Lou & Wang (2006) . 

In this subsection, we turn our attention to self-similar 
flows with shocks. From now on, subscripts 1 and 2 rep- 
resent upstream and downstream sides of a shock, respec- 
tively. In particular, we use Ci and C2 to represent Co of the 
upstream and downstream sides, respectively. Because it in- 
volves local sound speed with respect to the shock reference 
framework in both upstream and downstream sides, we also 
calculate the corresponding sound speed Cs = {'yP/pY^'^ = 
{'yCom^''^a^^^)^^^ for each branch of solutions. 

We begin with free-fall core collapse solutions. From the 
discussion of collapse solutions without crossing the sonic 
singular surface, any of this kind solution will cross the sonic 
singular surface even number of times, either smoothly or 
by shocks. By inspecting this topological characteristics and 
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considering the simplest case of a single shock, we infer that 
this type of solutions, with shock jumps across the sonic 
singular surface, should be also possible to cross the sonic 
critical line smoothly at some critical point. Based on this 
observation, we specify a type I eigensolution at a given 
sonic critical point and integrate away from it in both direc- 
tions. Let us take solutions shown in Figure [10] as examples 
of illustration. In the comoving reference framework of the 
shock, the outer part is supersonic and is thus the upstream 
side, and the inner part is the downstream side. Here we 
apply the matching procedure in the a — v phase diagram 
introduced by Hunter (1977). A notable differece from the 
case of 7 7^ 4/3 is that the value of Co and the shock posi- 
tion will affect the value of Co and thus the sonic singular 
surface in the other solution branch. We also have a con- 
siderable freedom to construct a shock in one solution at a 
chosen place and then integrate forward. Such a numerical 
solution may approach a certain asymptotic solution, or it 
may encounter the sonic singular surface. Remember that 
when a numerical integration is from the downstream side 
to the upstream side, the square of the Mach number M2 of 
the downstream side should be within the range of (0.125, 1) 
and thus the shock position Xs must be in a corresponding 
range. The Xa value in the example of Figure [10] is around 
0.4578. This kind of gravitational core collapse solutions to- 
gether with other solutions investigated previously, such as 
Shu (1977) and Lou & Shen (2004), may describe a possible 
stage of star formation in molecular clouds. 

Tsai & Hsu(1995), Shu et al. (2002) and Bian & 
Lou(2005) connected the outer singular isothermal sphere 
(SIS) solution with either LP type solution or free-fall so- 
lution in the inner region by shocks in an isothermal gas. 
Using the matching procedure in the a — v phase diagram, 
shock flow solution of this kind with the free-fall asymp- 
totic solution at small x also exists in our polytropic case 
of 7 = 4/3. This particular kind of shock solutions depicts 
the following scenario. Initially the outside gas is in a radial 
force balance and the collapse starts from the central core 
region. Effects such as changes in the centre propagates out- 
wards in the form of a self-similar shock. Materials are blown 
out by this shock. Because the gravity is stronger than the 
pressure force, materials eventually stop moving outwards 
and fall towards the centre. 

While all LP type asymptotic solutions degenerate to 
the Einstein-de Sitter solution in our case of 7 = 4/3, shock 
solutions can also be constructed to connect the inner Ein- 
stein de Sitter solution with an outer SPS (see Fig. I12|l . Nat- 
urally, the outer SPS part is the upstream side with v\ — 0, 
and the inner Einstein-de Sitter solution is the downstream 
side with V2 = 2x^/3 and 0:2 = 2(1 + 2^/3C2)'^/S at Xs 
where C2 is set to an appropriate value. Using 112, 02 and 
C2, we can express vi in terms of the scaling index a. The 
condition of vi — then appears as a quadratic equation 
of a, which needs to be solved for a with a < —2/3. Once 
this is done, we use V2, 0.2, C2 and the relevant root(s) of a 
to calculate the Mach number A^2 on the downstream side 
to check whether the requirement 1/8 < M2 < 1 is met. 
Once everything is complete and consistent, one of this kind 
of shock solutions is then constructed. We show an example 
here in Figure [T2l with C2 = 1.5, a = -0.799, 0L2 = 0.1252, 
leading to Ci = 0.383, ai = 0.0207 correspondingly, where 
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Figure 11. Examples of inner free fall and shock jumps to the 
SPS outer part and other asymptotic flow solutions far away. For 
the SPS shock connection, the shock is located at Xs = 1.23 with 
C\ = 0.6623 and C2 = 0.6696 and the inner solution crosses the 
sonic singular surface again at = 0.992 smoothly. The mass 
at the centre is r?i(0) = 1.3039 and the mass enclosed within the 
shock front is m{xs) = 1.3629. For the three outer asymptotic flow 
solutions from the top in order, we have relevant shock parameters 
{Ci, C2, Xs} to be {0.1810, 0.6690, 1.8}, {0.5020, 0.6690, 1.6}, 
{0.6276, 0.6690, 1.4}, respectively, and relevant flow parameters 
{H, L, K\ of asymptotic solution l|46] l to be 10.082. -0.1, -1.0}, 
{0.080, 0.0, -0.3}, {0.080, 0.0, -0.07}, respectively (see asymp- 
totic solution 146 l l. 
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Figure 12. This figure illustrates a special case (the heavy solid 
line) of shock solution which connects the inner Einstein-de Sit- 
ter solution and the outer SPS with 7 = 4/3. Here, C'2 = 1-58 
and the corresponding value of a is —0.80. This kind of solu- 
tions may be regarded as a limiting solution as C2 approaches 
1.58 while keeping a the same. Also shown in this Figure are 
different upstream solutions of C2 = 1-3, 1.4, 1.7, 1.8, respec- 
tively. For these four upstream solutions away from the SPS, 
the three relevant parameters {H, K, L} for asymptotic so- 
lution gg are {0.134, -0.007, 0.068}, {0.126, 0.00, 0.041}, 
{0.106, 0.00, -0.024} and {0.101, 0.00, -0.041}, respectively. 
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the self-similar shock position 2:3 can be any positive values 
because of the scale invariance property. 

Shocks can also be inserted to connect the only two an- 
alytic solutions available, namely, the static SPS solution 
outside and the Einstein-de Sitter solutions inside. To con- 
struct this kind of shock solutions, the upstream quantities 
are vi = and ai = ExV"^ with Ci (i.e., the upstream 
Co) satisfying conditions (|34|l and (|35|) for the existence of 
SPS, while the downstream quantities are V2 = 2xs/?> and 
02 = 2/[3(l + 2\/2>C2)] from the Einstein-de Sitter solution 
l|36p . With these constraints, it is straightforward to deter- 
mine C2 = 1.58, Ci = 0.42 and a = -0.80. Figure [12] shows 
this solution with the shock position at Xg = 2.0. It is easy 
to see from the figure that this solution represents the lim- 
iting solution of a solution family of C2 approaching 1.58 
with = 0. It is unlike the isothermal results of Tsai & 
Hsu (1995) where this kind of solutions is a limit of a family 
of breeze solutions (Shu et al. 2002). Instead, it is a critical 
state to distinguish asymtotic outflow and inflow solutions. 
For C2 being slightly larger than 1.58, the asymptotic solu- 
tion represents an inflow, while for C2 < 1.58 the asymptotic 
solution corresponds to an outflow. In fact, this Einstein-de 
Sitter shock model can be applied to an explosion process 
with a stellar interior as an alternative of the rebound shock 
model of Lou & Wang (2006, 2007) described at the begin- 
ning of the next paragraph. The major difference here is a 
constant density within the shock front instead of being a 
diverging density near the centre; outside the shock front, 
the density approaches a power-law scaling with either in- 
falling or outgoing stellar materials. When this shock front 
reaches the photosphere of the progenitor, we start to see 
observable effects of a supernova in optical bands. 

Lou & Wang (2006) utilized a self-similar polytropic 
model to construct the gravitational core collapse and re- 
bound shock processes in supernova explosions. Their con- 
ventional polytropic model solution with 7 7^ 4/3 is to con- 
nect quasi-static solutions at small x with outer asymptotic 
flow solutions at large x by outgoing shocks. That model 
was recently generalized to include a random magnetic field 
using a magnetohydrodynamic (MHD) approach and to ex- 
plore the origin of strong magnetic fields of compact ob- 
jects (Lou & Wang 2007; Wang & Lou 2007). Based on our 
model framework here, this can also be done for a general 
polytropic gas with 7 = 4/3 and thus q = 2/3. Starting nu- 
merical integrations both from the centre and from infinity 
(actually a sufficiently large x) and choosing a proper meet- 
ing point to match solutions in the a — v phase diagram, we 
adjust the shock position and parameters of outer asymp- 
totic solution (1461) to construct sensible solutions (e.g., Lou 
& Shen 2004). Alternatively, instead of the above matching 
procedure, we can also start a numerical integration from 
the vicinity of the centre and then choose a certain point as 
the shock location. Shock jump conditions (|71|l — H73|l de- 
termine all physical variables on the upstream side of a 
shock. A further numerical integration outwards until x is 
sufficiently large completes the solution construction proce- 
dure. Using the numerical solution thus obtained, we can 
match with asymptotic solutions to determine relevant pa- 
rameters. We find that the self-similar shock position can 
only exist within a finite interval of x (e.g., in the case of 
a — —0.68, the self-similar shock position falls within the 
range of 1.3 < Xs < 4.23). Outside this interval of x, the 
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Figure 13. Shock solutions with quasi-static asymptote of equa- 
tion II6II I and II62I I at small x are shown here. The power in- 
dex parameter a is —0.68. In the downstream side, we use the 
same quasi-static solution while two uptream branches are dif- 
ferent. The parameters in the downstream side for the quasi- 
static solution are C2 = 0.6623, B = 1.00, W = 0.10, V = 
0.06375, ^ = 1.8824. The upstream branch converges to asymp- 
totic solution 1 146 I I at large x. Thus each upstream solution 
can be identified by a set of parameters {Ci, Xs, H, L}, 
where Xs is the shock location and the last two are the co- 
efficients in asymptotic solution l|46| l. In this illustration the 
two sets of parameters are {0.62758, 2.5, 1.0096, -0.5323} and 
{0.62349, 2.7, 1.0149, -0.4486}. The dash-dotted curves are cor- 
responding segments of Vc curve defined by condition JSTJ. 



square of Mach number 7M2 on the downstream side would 
be smaller than 1/8 which is unphysical by our analysis on 
self-similar shock conditions. This 7 = 4/3 rebound shock 
model for supernovae may be more appropriate in certain 
aspects. During the initial phase for the emergence of a re- 
bound shock in the dense stellar core, neutrino pressure, ra- 
diation pressure and gas pressure together may be modelled 
by a polytropic mixture of 7 = 4/3. The diverging density 
near the centre is expected to create a highly degenerate 
core there. 

We also construct shocks connecting asymptotic solu- 
tions (|46|) at infinity and (|57[) near the centre. In the above 
analysis, we know that the inner part of this solution can 
only appear in the first quadrant in the plane of —v ver- 
sus X. The numerical treatment starts from the centre and 
goes outwards. Before the solution meets the sonic singular 
surface, jump conditions are included to introduce a shock 
across the sonic singular surface. We then continue to in- 
tegrate outwards until x is sufficiently large to match with 
asymptotic solutions at large x. Numerical experiments show 
that almost all such solutions match asymptotic solution 
(|46|) . in which mass density and flow velocity both converge. 
It is also possible that after a shock jumping across the sonic 
singular surface and integrating outwards, the radial flow 
velocity decreases rapidly so that it crashes onto the sonic 
singular surface again. In our numerical experiment, we do 
not find solutions with twin shocks or others, which jump 
across the singular surface twice or more (see Bian & Lou 
2005). 
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Figure 14. A similarity shock solution witii a = —0.68 and 
Ci = 0.85874, C2 = 0.9 is illustrated. It connects asymptotic 
solution 1146 l l with parameters H = 0.2585, L = —1.0301 and 
asymptotic solution II57I I with parameters R = —1.1235, = 
0.1, A = —2.9778. The solid curve represents the shock solution 
and the dashed curve represents corresponding segments of Vc 
curve defined by condition JSTJ. 



6 DISCUSSION AND CONCLUSIONS 

We have explored and examined the similarity flow solution 
structures in a general polytropic gas with a polytropic index 

7 = 4/3. Previously, Goldreich & Weber (1980) considered a 
special case of 7 = 4/3 with a constant specific entropy. By 
their assumptions and analysis, only homologous collapse 
solutions exist by invoking the time reversal invariance, i.e., 
the radial flow velocity u takes on the form of 2r/{3t) until 
the mass density vanishes at a certain point. Yahil (1983) 
mentioned the case of Goldreich & Weber (1980) as a spe- 
cial limit. In reference to earlier work and based on a self- 
similar transformation, we systematically examined the case 
of 7 = 4/3 with specific entropy conservation along stream- 
lines. We have substantially generalized the earlier analyses, 
discovered new asymptotic solutions, and constructed vari- 
ous self-similar solutions without or with shocks. 

In reference to earlier analyses of Goldreich & Weber 
(1980) and Yahil (1983), our model framework mainly fo- 
cuses on 7 — 4/3 with the conservation of specific entropy 
along streamlines, which is more general and perhaps, closer 
to reality than the conventional polytropic gas of a constant 
specific entropy everywhere at all times. Of course, the case 
of a constant specific entropy is also possible and can be 
properly accommodated and treated within our polytropic 
model framework of 7 = 4/3. Under our more general for- 
malism, we extend the work of Goldreich & Weber (1980) 
and obtain many interesting results. The solutions are di- 
vided into two broad classes: solutions with a = —2/3 pre- 
cisely belong to Class I and solutions with a < —2/3 belong 
to Class II. For the situation of —2/3 < a < as mentioned 
at the beginning of deducing asymptotic behaviours, the di- 
vergent velocity at large x is not of interest and hence we 
only consider two classes I and II solutions. 

Class I solutions are characterized by P cx p'^^^ with the 
proportional coefficient related to the specific entropy being 



an arbitrary function of x, while for Class II solutions, this 
proportional coefficient depends on the enclosed mass M in 
a power-law form. We discuss these two classes separately. 

6.1 Class I Self-Similar Solutions 

Class I self-similar solutions represent a substantial exten- 
sion of the special solutions with a constant entropy derived 
by Goldreich & Weber (1980). For an astrophysical system 
such as stars, the specific entropy is not expected to be a 
global constant in general. For a stellar interior, this de- 
pends on the competition between thermal kinetic energy 
and Fermi energy as determined by the mass density. Quali- 
tatively speaking, especially for a compact object, the closer 
to the centre, the closer the material is in a degenerate state; 
this would correspond to a smaller specific entropy. However, 
the density is relatively small and the temperature is rela- 
tively low in the outer part of a star, perhaps also leading 
to a lower level of entropy. We do not yet know the exact 
distribution of specific entropy within a star so far. Thus the 
case of a constant entropy is the simplest to consider and 
provides a certain sense for a homologous dynamic process. 
The model analysis of this paper is more general and allows 
for a fairly arbitrary distribution of specific entropy along 
streamlines. Meanwhile, the radial velocity profile remains 
always equal to 2r/(3t). For a given time t, the radial veloc- 
ity increases linearly with increasing radius r. Hence, this 
solution can be valid within a finite radial extent. It turns 
out that the mass density vanishes at some place referred as 
the outer boundary of the fiow system. 

According to the model analysis of Goldreich & Weber 
(1980), a pre-coUapse progenitor star of a static configura- 
tion may evolve into a homologous core collapsing phase (see 
Figure [T]) , when the pressure suddenly decreases by a frac- 
tion within a range of ~ 2.9%. Early simulations of Bethe et 
al. (1979) indicated a substantially larger pressure reduction 
of 26% is needed in order to initiate collapse in supernova 
explosions. The much smaller fraction change of pressure re- 
duction for a homologous core collapse given by Goldreich 
& Weber (1980) is actually related to the assumption of a 
constant specific entropy (i.e., their constant k) in space and 
time. Requiring specific entropy conservation along stream- 
lines and allowing the specific entropy to be a function of 
space and time, it is possible to have a homologous core 
collapse for a much larger fractional change of pressure re- 
duction. In our more general analysis and notations, we find 
that other forms of g{x) instead of g{x) = 1 can give rise 
to a fractional change of 26% or larger for a pressure reduc- 
tion. Physically, this corresponds to different distributions 
of specific entropy along streamlines. 

To illustrate this case specifically, we choose g{x) = 
1/(1 + ex) where e > is an adjustable parameter to gradu- 
ally modify the shape of g{x). When x is sufficiently small, 
g{x) is nearly equal to 1, analogous to the g{x) = 1 case. 
Globally g{x) is a decreasing function with increasing x. 
When e is small, g{x) decreases slowly and only deviates 
from g{x) — 1 case when x is sufficiently large. For large 
e > 0, the result will differ considerably from that of Gol- 
dreich & Weber (1980). We carry out such a g{x) experiment 
numerically. 

Substituting dimensional quantities into the dimension- 
less state function p = g{x)a'''^''^ , the dimensional equation 
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Table 2. When g{x) takes the form of 1/(1 + ex), the results 
are shown in Figure 1151 As /(O) is sufficiently large, the outer 
boundary at a finite x where f{x) = is extremely small and 
thus g{x) can be almost treated as a constant; these results vary 
little as compared with those of Goldreich & Weber (1980) for 
/(O) — > +00 corresponding to the Lane-Emden equation. Results 
in Figure [TSl mainly focus on the limiting case of /(O) f^. For 
possible f{x) solutions of a homologous core collapse, we use Tp 
to denote the range for fractional change of pressure variation as 
/(O) increases from fc to infinity. 



e 


fc 


p/pc at fc 


rp 





4.67047 


0.00655 


2.9% 


0.01 


4.58642 


0.00693 


3.6% 


0.05 


4.28755 


0.00848 


6.4% 


0.1 


3.98354 


0.0106 


10% 


0.3 


3.23130 


0.0198 


26% 


0.5 


2.83696 


0.0293 


42% 




Figure 15. Results of the normalized f{x)/f(0) is displayed for 
g{x) = 1/(1 + O.lx) and e = 0.1. In this case, the minimum of 
/(O) is fc = 3.9835 for physical solutions. 

of state can be written explicitly as 

P=£M(4^G)i/V/' , (84) 
yi 

where for a given time t, g(x) corresponds to a radial dis- 
tribution of specific entropy. Parameter A also varies for 
different values of chosen /(O) for a certain system in which 
the total enclosed mass M is conserved and thus the value of 
A parameter can be deduced from equation (|23p . The vari- 
ation in A value actually corresponds to the variation range 
for fractional change in pressure denoted by rp \ in Goldreich 
& Weber (1980), this Tp is 2.9% as /(O) increases from fc to 
infinity. For any given form of g(x), the limiting case is the 
Lane-Emden equation (e.g., Chandrasekhar 1939) as long as 
g{x) ^ 1 as a; 0^. Of course, this applies to our chosen 
form of g{x) = 1/(1 + ex) as a; — > 0"*". 

Tableland Figure [15] show major results for a range of 
e > values. Numerical experiments indicate that as e in- 
creases from 0, the limiting value fc of /(O) decreases while 
the density ratio p/pc increases. More importantly, the range 
of fractional change rp by which the pressure can be reduced 



for a homologous core collapse becomes larger and larger. 
Table [2] shows that for e = 0.3, we have rp = 26%; in other 
words, for such a large fractional change in pressure reduc- 
tion (e.g., Bethe et al. 1979) in order to initiate supernovae, 
it is still possible for a homologous core collapse prior to 
the development of a rebound shock. It is conceptually im- 
portant that a different specific entropy distribution of g{x) 
from a constant value can lead to a better agreement with 
numerical simulations; this appears more effective than the 
inclusion of a less massive core in the centre as mentioned 
in Goldreich & Weber(1980). 

6.2 Class II Self-Similar Solutions 

In addition to extensions of Goldreich & Weber (1980) dis- 
cussed in the above subsection, we also substantially gener- 
alize the self-similar solution space for 7 = 4/3 by adjusting 
the scaling index a. In contrast to 7 7^ 4/3, a straightforward 
analysis with 7 = 4/3 leads to an exact value of g = 2/3 that 
is independent of a. The dimensional equation of state then 
takes the form of P oc M^^'^p'*^^ with a constant propor- 
tional coefficient. We can compare the thermal energy ksT, 
where is the Boltzmann constant, with the Fermi energy 
ef- Neglecting the rest mass in the relativistic regime, the 
relationship between total energy e and momentum p for 
a single particle can be written as e = cp where c is the 
speed of light, leading to Ef oc p^^^ . In our model, the state 
function gives ksT oc P/p oc M^/^p^^^. It follows that 

^ oc . (85) 

Ef 

The enclosed mass M is a non-decreasing function in radius 
r. At a given time t, we see from this relation that, at small 
r, the enclosed mass is small and hence this ratio is also 
small. Physically in the inner core of a star, where materials 
are highly condensed and may be close to a degenerate state, 
the specific entropy is low. 

A modified self-similar transformation is introduced for 
7 = 4/3. In the self-similar transformation for a conven- 
tional polytropic gas (i.e., P = Acp^ with globally constant 
K, at all times), the sound speed appears either explicitly or 
implicitly. In contrast, we here use an integration constant 
Co relating to the sound speed and transformation ((Tjl does 
not involve the sound speed because of the uniqueness of 
the 7 = 4/3 case; this Co coefficient is in fact allowed by 
the transformation and is an adjustable parameter in our 
analysis for astrophysical applications. At a deeper level, we 
realize that as the special self-similar transformation ([7]) for 
7 = 4/3 does not involve the sound speed, we have a scal- 
ing invariance (|33|) which simplifies our theoretical analysis 
considerably. 

By comparisons and analogies of solutions known for 
7 7^ 4/3, we try our best to find the counterpart solu- 
tions and to discover new solutions for 7 = 4/3. Global 
analytic solutions, i.e., static SPS solution (|34p and (|35|l 
and Einstein-de Sitter expansion solution (|36|l , still exist for 
7 = 4/3 with some modifications. Analytic asymptotic so- 
lutions of various kinds are also derived for both large and 
small X. However, the LP- type solution no longer exists for 
7 = 4/3 (thus q = 2/3) except for rare situations, while other 
counterpart solutions are readily found. In addition, a new 
type asymptotic solution (|57|) is discovered in the regime of 
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small X. It seems that this type of asymptotic polytropic 
solutions only exists for 7 = 4/3. 

We have also examined properties of the sonic singu- 
lar surface and the sonic critical line of coupled nonlinear 
ODEs (|3H) . A salient feature of 7 = 4/3 case is that all 
sonic critical lines are straight and pass through the origin 
a:: = and v = m the —v(x) versus x presentation; while 
first revealed by extensive numerical experiments, these re- 
markable results can be proven analytically. It is also fairly 
straightforward to derive two eigensolutions to smoothly 
cross the straight sonic critical line. In later analyses, we 
realize that the sonic critical line cannot fully represents the 
behaviours of the sonic singular surface, especially for con- 
structing shocks. We hence use Vc{x) defined by equation 
(|8ip for each solution, which is another curve on the sonic 
singular surface and tightly relates to the current solution, 
to show the interrelation between the current solution and 
the sonic singular surface. According to definition (|8ip . the 
solution meets the singular surface if and only if the solution 
and the corresponding Vc intersects; meanwhile, a shock so- 
lution jumping across the sonic singular surface also jumps 
across the corresponding Vc in v{x) versus x plane. The stan- 
dard Runge-Kutta scheme (e.g.. Press et al. 1986) is used to 
numerically integrate coupled nonlinear ODEs (|31|l and (|32p 
to connect various asymptotic solutions and eigensolutions 
across the sonic critical line. We also construct possible so- 
lutions with and without shocks for potential astrophysical 
applications. From the behaviours of these semi-complete 
solutions, we can see that many such solutions have similar 
behaviours as those for 7 7^ 4/3 in a qualitative manner. And 
our solutions can be sensibly regarded as limits of 7 —> 4/3 
when a polytropic gas becomes relativistically hot or degen- 
erate. Our analysis for 7 = 4/3 does share a certain common 
characteristics with cases for 7 7^ 4/3. 

We can readily solve for eigensolutions crossing the 
straight sonic critical line. Once the slope k of the sonic 
critical line is positive, we can construct a new type of 
self-similar solution characterized by an expanding central 
spherical void within which the enclosed mass is zero or neg- 
ligible. Around the edge of such an expanding void, there ex- 
ists an overdensed shell where the density variation becomes 
rather steep. Diffusion processes are expected to smooth out 
such relatively steep gradients locally. To consider proper- 
ties of spherical void boundary, we first note that such a void 
expands with a radial speed 

re = —ar^jt re oc f"" , (86) 

where re stands for the radius of the spherical void edge. The 
spherical void edge evolves as a power law of time t with a 
scaling index —a. One also notes by equation (|82|) that the 
density gradient approaches a negative infinity near the void 
edge. We expect that in a narrow region near the spherical 
void edge, materials are actually diffused instead of being 
so sharply distributed as shown by our solution mathemati- 
cally. Within this narrow region, the local evolution does not 
behave self-similarly and may not be spherically symmetric, 
while the overall self-similar profile remains on large scales. 
In the outer part, the mass density scales as p oc r'^''" and 
the radial flow velocity remains finite with a wind. In fact, 
it is also possible to construct various shock solutions with 
a central void. 

At this stage, we may outline a physical scenario in the 



context of a supernova explosion. During the core collapse 
of a progenitor, neutrons are formed in abundance and neu- 
trinos of relativistic energies are released. In a high-density 
environment, neutrino opacity is extremely high so that neu- 
trino pressure, radiation pressure and gas pressure work to- 
gether to drive the central core expansion. In the relativistic 
regime, we may ignore tiny neutrino masses and regard the 
neutrino gas as polytropic with an index 7„ = 4/3. Simi- 
larly, the radiation pressure resulting from the photon gas 
trapped in the stellar interior can also be regarded as poly- 
tropic with an index 7^ = 4/3. In the hot stellar core of high 
temperatures, we may approximate the thermal gas pressure 
as polytropic with an index 7p = 4/3. It might be conceiv- 
able that under certain situations, the neutrino pressure is 
so overwhelming such that a central void may start to form. 
As the outer part expands and density drops, neutrinos es- 
cape while the radiation and thermal gas pressures continue 
to drive the expansion. It should be emphasized that in real 
situations, a grossly spherical void may still encompass ma- 
terials here and there but the mass density inside is substan- 
tially lower than that of surroundings. 

In this context, we note the model work of Fillmore & 
Goldreich (1984b) who considered a collection of collision- 
less particles in an expanding universe of the Einstein-de Sit- 
ter form. There is no pressure effect from particles in their 
model. In essence, the background Einstein-de Sitter expan- 
sion prescribed is similar to the rapid expansion driven by 
the thermal pressure force in our model, both providing the 
tendency for particles to move outwards in competition with 
the inward self-gravity. The key physical difference is that 
the Einstein-de Sitter expansion of the universe is homoge- 
neous (presumably driven by the ubiquitous dark energy) 
while our gas expansion is driven by the thermal gas pres- 
sure closely related to gas mass density and temperature. 
Not only in the case 7 = 4/3, self-similar void solutions can 
also be constructed for 7 7^ 4/3 and isothermal gas which 
we shall investigate more thoroughly in separate papers. 

Besides certain similarities with previous polytropic 
model analysis with 7 7^ 4/3, the case of 7 = 4/3 carries 
its own unique features. First, because of scaling invariance 
H33|) . various self-similar solutions can be readily classified, 
especially for the two eigensolutions across the sonic criti- 
cal line. Once solution properties at a chosen point x have 
been examined completely, other points will have the same 
solution characteristics by scaling invariance (I33p . This sim- 
plifies the analysis to a considerable extent. Fundamentally, 
the cause of this scale invariance (|33} is due to the fact that 
the sound speed is not involved in self-similar transformation 
0. In various solutions, the case of 7 = 4/3 also shows some 
differences: (i) the LP type solution does not exist, except 
for rare situations (see Appendix[B} ; (ii) when discussing the 
quasi-static solution at small x, two sensible roots may be 
found for 7 < 4/3 (see Lou & Wang 2006), while one of the 
two roots is always unacceptable for 7 = 4/3 with only one 
sensible root being available in our model calculations; (iii) 
it is no longer possible for a quasi-static solution at small 
X to show a vibration behaviour here. The solution quickly 
converges to an outer asymptotic solution; and (iv) the sonic 
critical lines with constant density are straight lines emanat- 
ing from the origin in the —v(x) versus x presentation. 
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6.3 Conclusions 

To sum up this paper, we iiave explored possible self-similar 
solutions, both analytical and numerical, for a generalized 
polytropic gas with 7 — 4/3. The classical analysis of Gol- 
dreich & Weber (1980) for a constant specific entropy ev- 
erywhere at all times is substantially extended by specific 
entropy conservation along streamlines with specific entropy 
dependent on time and space. This difi'ers from what Yahil 
(1983) did in this context. In addition to counterparts of 
various previously known types of polytropic solutions with 
7 7^ 4/3, we find two new asymptotic solutions. One notable 
feature is that sonic critical lines are all straight lines em- 
anating from the origin in the —v{x) versus x presentation 
with constant densities. Using all asymptotic solutions avail- 
able, two eigensolutions across the sonic critical line and self- 
similar shock conditions, global semi-complete solutions are 
constructed numerically. Two classes of self-similar solutions 
are investigated separately according to the value of the scal- 
ing index a. For Class I solutions with a — —2/3 precisely, 
we can simulate the homologous evolution of a fiow system 
once the distribution of specific entropy is prescribed. These 
more general solutions for a homologous core collapse (Gol- 
dreich & Weber 1980) may be utilized to model the dynamic 
formation of an inner compact core from a pre-coUapse stel- 
lar interior. Collapsing solutions of Class II with a < —2/3 
may also explain the formation of compact objects and other 
similar fiow systems, while expansion solutions with shocks 
can be utilized to model supernova explosions (e.g., Lou & 
Wang 2006, 2007). 

By equations (|8]) — (|12|) , we emphasize several aspects of 
solutions for (3a -I- 2) = and for (3a -I- 2) 0. Based on our 
analysis, the Einstein-de Sitter solution with 7 = 4/3 exists 
for (3o 2) 7^ 0, (3a + 2) ^ and (3a + 2) = 0. Except 
for this special Einstein-de Sitter solution. Class I solutions 
valid for (3a + 2) = cannot be obtained by taking the 
limit of (3a -I- 2) — > for Class II solutions. In other words. 
Class I and II solutions are qualitatively different solutions 
and we need to consider them separately. By equations ([8]) 
and (|25[) during the limiting process of (3a -I- 2) ^ 0, we 
must require {ax + «) ^ in order to have a finite reduced 
mass m{x). Only the Einstein-de Sitter solution and Class I 
solutions with (3a -I- 2) = bear this unique feature for the 
reduced fiow speed v{x) while all other Class II solutions are 
excluded by this limiting procedure. 

In the course of investigation, we realize the possibil- 
ity of constructing self-similar solutions for dynamic evolu- 
tion of central spherical void in a fiow system involving self- 
gravity and thermal pressure. Here, the thermal pressure 
force drives the gas expansion sufficiently fast and creates 
a central spherical void by pushing materials outwards. By 
specific examples, we now prove by analytical and numerical 
calculations that a spherical void can indeed form in astro- 
physical fiow systems under the joint action of thermal pres- 
sure force and self-gravity. We expect that such processes 
could happen in association with supernova explosions and 
evolution of supernova remnants. 

At the beginning of our model formulation, several 
physical effects, such as nuclear reactions, radiation pres- 
sure, neutrino transport, general relativistic effects, rota- 
tional effects and magnetic field, are not taken into account. 
Under various situations, these effects can be very impor- 



tant in real astrophysical systems. Given these approxima- 
tions and idealizations of our model, it is still hoped that 
this simple theoretical model framework may catch certain 
essential characteristics or features of fiow phenomena of rel- 
evant scenarios and interpretations. 
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(|12|) is equivaleni io equations (|8} and (|9]). Assuming Tay- 
lor series expansions in ihe vicinity of s = 0, we write the 
solutions as 



(Bl) 



where Vk and a; are constant coefficients with cto 7^ 0. By 
equation (|9}, the enclosed mass is given by 



i{x) = I y^a{y)dy = ^ 

"'0 ,_n 



(; + 3)- 



(B2) 



While from equations and ([U}, we have 

[ax + v) 2 



(3a + 2) 



The two expressions of m{x) should be equal, giving rise to 
a series of relations among the coefficients of a{x) and v{x), 



vo^O , vi^ 2/3 
4(a -I- vi)ai + 4u2ao = (3a -|- 2)qi 



(B4) 
(B5) 



APPENDIX A: THE EXPLICIT FORM OF V 

AND a' 

Using Cramer's rule in equations (|31|l and (I32p . one can 
easily deduce the explicit forms of v'{x) and a'(x), namely, 



v'ix) = Vix)/V{x) , 
a{x)/a{x) = A{x)/V{x) , 

where 

N (ax + v)'^ , . 
V{x) = -^^^5 — -wrCi + {a+ l)(ax + v)v 



(Al) 
(A2) 



(3a + 2) 

f 3a- 



A{x) = 2{ax + v) (l - ^) 



V \ {ax + v) 
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a — (a + l)v 
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(3a + 2) 

2Co / ax + v ^ ^4/3^ 
3 V3a + 2 



D(a;) = {ax + v) 



2 4Co 4/3 (ax + v 



V3a + 2 



2/3 



(A4) 



(A5) 



The sonic singular surface corresponds to ©(a;) = 0; together 
with either V{x) = or ^(a;) = 0, the sonic critical line 
is then determined. Note that V(x) = is equivalent to 
A{x) = on the sonic singular surface. 



APPENDIX B: EXISTENCE CONDITION FOR 
LARSON-PENSTON TYPE SOLUTIONS 

Starting from equations (|8|l — (1121 ). we briefly discuss the ex- 
istence of Larson-Penston (LP) type solutions for a general 
case of a 7^ —2/3 without constraining 7. In fact, equation 



Besides, the reduced pressure can also be written as a series 
expansion in the form of 



Com'a^ — Cqx 



3q 



E 



-X 



{i + 3y 



E 



CtkX 



(B6) 



Substituting all these series expansions into equation (|10[) 
and comparing coefficients of the same powers of 2;, we have 
the following conclusions. For an arbitrary q in general, con- 
sider the power factor x^'' in equation (IB6p . If 3q is not an 
integer, the power index of every term in p is not an integer, 
and the terms thus cannot have the same power of other 
terms in equation pOjl . Consequently, no such a series so- 
lution exists. A necessary condition for the existence of LP 
type solution is that q takes the form of J/3 with J being 
an integer. For example, J = 0, q = and thus a = 7 — 2 
and Co = 1, we can readily obtain the following asymp- 
totic solution (equations 28a and 28b in Suto & Silk 1988), 
namely 
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ao-- 



a+'^]x^ + - 



a{x) 



Qo 



67 



Qo 



(B7) 



(B8) 



For J = 2, g = 2/3 and p = Com'^^^a*^^ , one readily 
obtains 



uo 



: 



vi = 2/3 



V2 







(B9) 
(BIO) 



Qo = I (l + 2v^Co) \ Qi = . 

For any given index integer fc > 1 in the series expansion, if 
we have already determined coefficients ai and Vi+i where 
0<i<k — 1, a comparison of the coefficients of each side 
of ODEs (H} and pop will give a pair of linear equations for 
Qfc and Vk+i, which has a unique solution for Ok and v^+i- 
Thus, coefficients q; and Vi {i > 1) have only one solution. 
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On the other hand, Vi = and ai = (i > 1) gives a 
solution to this problem. Consequently for q — 2/3, only the 
Einstein-de Sitter solution exists and no LP type of solution 
can be found. One possible yet rare exception occurs when 
the coefficient determinant of one of the linear equations for 
Vi and ai becomes zero. We do not give more calculations 
on these special cases in this paper. 



APPENDIX C: REQUIREMENT ON Co FOR 
THE EXISTENCE OF ASYMPTOTIC 
SOLUTION (59) 

We denote the left-hand side of equation (I59|l by h{R) where 
f? < 2/3 is required by A < 0. For a sufficiently large value 
of |_R| with R<0,we have h{R) > 0; and for R 2/3", we 
also have a positive h{R). Taking the first derivative of h{R) 
and setting it equal to 0, we obtain only one root denoted 
by Ro, namely 

Ro = {6Cof^^{Sa + 2)-a (CI) 

for the minimum of h{R). Therefore if h{Ro) is also larger 
than 0, then equation h{R) = has no real roots and hence 
asymptotic solution (|59|l does not exist. On the other hand, 
if h{Ro) is smaller than zero, then equation h{R) — al- 
ways has two real roots. The critical case of h{Ro) = 
corresponds to a double root, and a critical value of Co = 
2*/^/6 « 0.4200 is thus known. 



